file <- paste0(ref_path, "/PFIRS/PFIRS 2019-2020 CT pulled 2022.xlsx")
xy_old <- read_excel(file)xy_old %>% head()## # A tibble: 6 × 18
## `Burn Date` `Burn Unit` `Acres Burned` `Total Tons`
## <dttm> <chr> <dbl> <chr>
## 1 2019-03-24 00:00:00 Office 0.01 5.4300000000000001E-2
## 2 2019-12-16 00:00:00 Mammoth Bar 0.01 0.1135
## 3 2019-04-16 00:00:00 Foresthill Loop 0.02 9.6000000000000002E-2
## 4 2019-02-08 00:00:00 Mammoth Bar 0.05 0.5675
## 5 2019-03-13 00:00:00 Foresthill Loop 0.05 0.24
## 6 2019-03-27 00:00:00 Foresthill Loop 0.05 0.24
## # ℹ 14 more variables: `Estimated Emissions` <chr>, Agency <chr>,
## # `Agency Unit` <chr>, Landowner <chr>, `Fuel Type` <chr>, `Burn Type` <chr>,
## # Latitude <chr>, Longitude <chr>, `Legal Location` <chr>,
## # `Min Elevation` <dbl>, `Max Elevation` <dbl>, County <chr>,
## # `Air District` <chr>, `Air Basin` <chr>xy_old <- xy_old %>%
rename(Burn_Date = `Burn Date`) %>%
rename(Burn_Unit = `Burn Unit`) %>%
rename(Acres_Burned = `Acres Burned`) %>%
rename(Fuel_Type = `Fuel Type`) %>%
rename(Total_Tons = `Total Tons`) %>%
rename(Burn_Type = `Burn Type`) xy_old <- xy_old %>%
filter(!Longitude == "UNK") %>%
mutate(Longitude = as.numeric(Longitude)) %>%
mutate(Latitude = as.numeric(Latitude))xy_old <- xy_old %>%
mutate(Longitude = ifelse(Longitude > 0, Longitude*(-1), Longitude))xy_old <- xy_old %>%
mutate(Year = year(Burn_Date))xy_old_2019 <- xy_old %>%
filter(Year == 2019)xy_old_2019 %>% head()## # A tibble: 6 × 19
## Burn_Date Burn_Unit Acres_Burned Total_Tons `Estimated Emissions`
## <dttm> <chr> <dbl> <chr> <chr>
## 1 2019-03-24 00:00:00 Office 0.01 5.4300000… 0.25
## 2 2019-12-16 00:00:00 Mammoth Bar 0.01 0.1135 0.32
## 3 2019-04-16 00:00:00 Foresthill … 0.02 9.6000000… 0.6
## 4 2019-02-08 00:00:00 Mammoth Bar 0.05 0.5675 0.32
## 5 2019-03-13 00:00:00 Foresthill … 0.05 0.24 0.6
## 6 2019-03-27 00:00:00 Foresthill … 0.05 0.24 0.6
## # ℹ 14 more variables: Agency <chr>, `Agency Unit` <chr>, Landowner <chr>,
## # Fuel_Type <chr>, Burn_Type <chr>, Latitude <dbl>, Longitude <dbl>,
## # `Legal Location` <chr>, `Min Elevation` <dbl>, `Max Elevation` <dbl>,
## # County <chr>, `Air District` <chr>, `Air Basin` <chr>, Year <dbl>From Jason Branz: “Typically if there are multiple records with the same date, burn name, and acres, it’s a duplicate.”
xy_old_2019 %>%
select(Burn_Date, Burn_Unit, Acres_Burned) %>%
filter(duplicated(.))## # A tibble: 112 × 3
## Burn_Date Burn_Unit Acres_Burned
## <dttm> <chr> <dbl>
## 1 2019-12-05 00:00:00 SPRA 0.25
## 2 2019-12-10 00:00:00 Unit 2 0.25
## 3 2019-12-10 00:00:00 SPRA 0.25
## 4 2019-12-11 00:00:00 Unit 4 0.25
## 5 2019-12-16 00:00:00 Big Chico Creek Ecological Reserve 0.25
## 6 2019-12-16 00:00:00 Big Chico Creek Ecological Reserve 0.25
## 7 2019-12-17 00:00:00 BCCER Prop Management 0.25
## 8 2019-12-17 00:00:00 BCCER Prop Management 0.25
## 9 2019-01-29 00:00:00 SPRA 1
## 10 2019-01-29 00:00:00 SPRA 1
## # ℹ 102 more rowsxy_old_2019 <- xy_old_2019 %>%
distinct(Burn_Date, Burn_Unit, Acres_Burned, .keep_all = T)sf_old_2019 <- st_as_sf(xy_old_2019, coords = c("Longitude", "Latitude"), crs = 4326)st_write(sf_old_2019, dsn = "shapefiles", layer = "PFIRS_2019_pulled2022", driver = "ESRI Shapefile", delete_layer = T)## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver## Deleting layer `PFIRS_2019_pulled2022' using driver `ESRI Shapefile'
## Writing layer `PFIRS_2019_pulled2022' to data source
## `shapefiles' using driver `ESRI Shapefile'## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 6: Field Burn_Dt create as date
## field, though DateTime requested.## Writing 1100 features with 17 fields and geometry type Point.## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: One or several characters
## couldn't be converted correctly from UTF-8 to ISO-8859-1. This warning will
## not be emitted anymore.xy_new <- read_excel(paste0(ref_path, "/PFIRS/PFIRS 2017-2022 pulled 2023.xlsx"))## Warning: Expecting date in A8155 / R8155C1: got 'NA'xy_new## # A tibble: 10,650 × 9
## `Burn Date` `Burn Unit` Agency `Acres Burned` Latitude Longitude
## <dttm> <chr> <chr> <dbl> <dbl> <dbl>
## 1 2017-01-02 00:00:00 PP Ridgetop Unit Cal F… 10 38.8 -123.
## 2 2017-01-03 00:00:00 Fairway/Bunker Calif… 2.5 39.2 -120.
## 3 2017-01-03 00:00:00 PP Ridgetop Unit Cal F… 10 38.8 -123.
## 4 2017-01-03 00:00:00 FIA Forest hill Priva… 10 39.0 -120.
## 5 2017-01-03 00:00:00 28A US Fo… 51 37.7 -119.
## 6 2017-01-04 00:00:00 Plum 52 US Fo… 1 39.5 -121.
## 7 2017-01-04 00:00:00 A-01 US Fo… 5 37.8 -119.
## 8 2017-01-04 00:00:00 PP Ridgetop Unit Cal F… 10 38.8 -123.
## 9 2017-01-04 00:00:00 5 US Fo… 59 37.6 -119.
## 10 2017-01-05 00:00:00 Big Chico Creek… CSU C… 0.1 39.8 -122.
## # ℹ 10,640 more rows
## # ℹ 3 more variables: `Fuel Type` <chr>, `Total Tons` <chr>, `Burn Type` <chr>xy_new <- xy_new %>%
rename(Burn_Date = `Burn Date`) %>%
rename(Burn_Unit = `Burn Unit`) %>%
rename(Acres_Burned = `Acres Burned`) %>%
rename(Fuel_Type = `Fuel Type`) %>%
rename(Total_Tons = `Total Tons`) %>%
rename(Burn_Type = `Burn Type`) xy_new <- xy_new %>%
mutate(Longitude = ifelse(Longitude > 0, Longitude*(-1), Longitude))xy_new <- xy_new %>%
mutate(Year = year(Burn_Date))xy_new <- xy_new %>%
distinct(Burn_Date, Burn_Unit, Acres_Burned, .keep_all = T)xy_new_2019 <- xy_new %>%
filter(Year == 2019)xy_new_2019 %>% head()## # A tibble: 6 × 10
## Burn_Date Burn_Unit Agency Acres_Burned Latitude Longitude Fuel_Type
## <dttm> <chr> <chr> <dbl> <dbl> <dbl> <chr>
## 1 2019-01-01 00:00:00 Forest M… Leoni… 1 38.6 -121. Slash
## 2 2019-01-01 00:00:00 Leoni Me… Leoni… 2 38.6 -121. Natural
## 3 2019-01-02 00:00:00 Foresthi… Calif… 0.1 39.0 -121. Brush
## 4 2019-01-02 00:00:00 Forest M… Leoni… 1 38.6 -121. Slash
## 5 2019-01-02 00:00:00 AIRPERSON Sierr… 2 38.9 -121. Slash
## 6 2019-01-02 00:00:00 Wyntoon … Wynto… 2 41.2 -122. UNK
## # ℹ 3 more variables: Total_Tons <chr>, Burn_Type <chr>, Year <dbl>sf_new_2019 <- st_as_sf(xy_new_2019, coords = c("Longitude", "Latitude"), crs = 4326)st_write(sf_new_2019, "shapefiles", "PFIRS_2019_pulled2023", driver = "ESRI Shapefile", delete_layer = T)## Deleting layer `PFIRS_2019_pulled2023' using driver `ESRI Shapefile'
## Writing layer `PFIRS_2019_pulled2023' to data source
## `shapefiles' using driver `ESRI Shapefile'
## Writing 1288 features with 8 fields and geometry type Point.sf_new <- st_as_sf(xy_new, coords = c("Longitude", "Latitude"), crs = 4326)st_write(sf_new, "shapefiles", "PFIRS_2017-2022_pulled2023", driver = "ESRI Shapefile", delete_layer = T)## Deleting layer `PFIRS_2017-2022_pulled2023' using driver `ESRI Shapefile'
## Writing layer `PFIRS_2017-2022_pulled2023' to data source
## `shapefiles' using driver `ESRI Shapefile'
## Writing 10252 features with 8 fields and geometry type Point.save(sf_new, file = "PFIRS_2017-2022_pulled2023.Rdata")CA <- st_read(dsn = paste0(ref_path, "/CA boundary/ca-state-boundary/", layer = "CA_State_TIGER2016.shp"))## Reading layer `CA_State_TIGER2016' from data source
## `C:\Users\ctubbesi\OneDrive - California Air Resources Board\Documents\Reference data\CA boundary\ca-state-boundary\CA_State_TIGER2016.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13857270 ymin: 3832931 xmax: -12705030 ymax: 5162404
## Projected CRS: WGS 84 / Pseudo-MercatorCA <- st_transform(CA, st_crs(sf_new))sf_new <- st_intersection(sf_new, CA)## Warning in st_is_longlat(x): bounding box has potentially an invalid value
## range for longlat data
## Warning in st_is_longlat(x): bounding box has potentially an invalid value
## range for longlat data## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesload(file = "Rdata/NF.Rdata")NF <- st_transform(NF, st_crs(sf_new))sf_use_s2(FALSE)## Spherical geometry (s2) switched offsf_new_noNF <- st_difference(sf_new, NF)## although coordinates are longitude/latitude, st_difference assumes that they
## are planar## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesmapview::mapview(list(NF, sf_new_noNF), col.regions=list("red","blue"),col=list("red","blue"))